Big Data and Machine Learning en el Mercado Inmobiliario
Clase 05: Prepocesamiento de un conjunto de datos (II)
Eduard Fernando Martínez González
Universidad de los Andes | BD-ML[1.] Imágenes satelitales
1.1 Motivación
Night lights and economic growth
Asian Financial Crisis: Java, Indonesia
Tomado de: Measuring Economic Growth from Outer Space
Genocide Event: Rwanda
Tomado de: Measuring Economic Growth from Outer Space
Puede acceder al documento aquí: Measuring Economic Growth from Outer Space
Night lights and Covid-19
Puede acceder al documento aquí: COVID19 y actividad económica: demanda de energía y luminosidad ante la emergencia
Puede acceder al documento aquí: Examining the economic impact of COVID-19 in India through daily electricity consumption and nighttime light intensity
Measuring the size and growth of cities using nighttime light
Puede acceder al documento aquí: Measuring the size and growth of cities using nighttime light
Housing Vacancy Rate (VHR) y Luces Nocturnas
Puede acceder al documento aquí: An estimation of housing vacancy rate using NPP-VIIRS night-time light data and OpenStreetMap data
CAUSES OF SPRAWL: A PORTRAIT FROM SPACE
Urban Land in Miami, FL
Tomado de: CAUSES OF SPRAWL: A PORTRAIT FROM SPACE
Puede acceder al documento aquí: CAUSES OF SPRAWL: A PORTRAIT FROM SPACE
1.2 ¿Datos raster?
¿Qué es un raster?
Tomado de: https://www.neonscience.org
Resolución
Tomado de: https://www.neonscience.org
Bandas de un raster
Tomado de: https://www.neonscience.org
1.3 Fuentes:
1.3.1 Night lights
DMSPL (1992-2013)
Puede acceder a los datos anuales aquí 1992-2103
VIIRS (2012-Currently)
Puede acceder a la web del proyecto aquí: link
Puede acceder a los datos mensuales aquí: link
Puede acceder a los datos diarios aquí: link
Harmonized (1992-2020)
Puede acceder a los datos y al paper aquí: link
1.3.2 Land cover
1.3.3 NASA
Puede acceder a datos de lluvias, temperatura, contaminación, vientos y otros en la página oficial del proyecto GES DISC de la NASA aquí: https://disc.gsfc.nasa.gov.
[2.] Manipular datos raster en R
0.0 Configuración inicial
Para replicar las secciones de esta clase, primero debe descargar el
siguiente proyecto
de R y abrir el archivo clase-05.Rproj.
Instalar/llamar las librerías de la clase
## Instalar/llamar las librerías de la clase
require(pacman)
p_load(tidyverse,rio,skimr,viridis,osmdata,
raster,stars, ## datos raster
ggmap, ## get_stamenmap
spatstat, ## analis de clusters
sf, ## Leer/escribir/manipular datos espaciales
leaflet) ## Visualizaciones dinámicas
## solucionar conflictos de funciones
select = dplyr::select2.1 Leer un raster
methods(class = "stars")## [1] [ [[<- [<- %in%
## [5] $<- adrop aggregate aperm
## [9] as_tibble as.data.frame as.owin c
## [13] coerce contour cut dim
## [17] dimnames dimnames<- droplevels filter
## [21] hist image initialize is.na
## [25] Math merge mutate Ops
## [29] plot predict print pull
## [33] replace_na select show slice
## [37] slotsFromS3 split st_apply st_area
## [41] st_as_sf st_as_sfc st_as_stars st_bbox
## [45] st_coordinates st_crop st_crs st_crs<-
## [49] st_dimensions st_dimensions<- st_extract st_geometry
## [53] st_interpolate_aw st_intersects st_join st_mosaic
## [57] st_normalize st_redimension st_sample st_set_bbox
## [61] st_transform_proj st_transform transmute write_stars
## see '?methods' for accessing help and source code## importar raster de luces: raster
luces_r = raster('input/raster/night_light_202301.tif')
luces_r## class : RasterLayer
## dimensions : 2990, 2919, 8727810 (nrow, ncol, ncell)
## resolution : 0.004166667, 0.004166667 (x, y)
## extent : -79.01042, -66.84792, 0.002082733, 12.46042 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## source : night_light_202301.tif
## names : night_light_202301
## values : -0.64, 2208.49 (min, max)## importar raster de luces: stars
luces_s = read_stars("input/raster/night_light_202301.tif")
luces_s## stars object with 2 dimensions and 1 attribute
## attribute(s), summary of first 1e+05 cells:
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## night_light_202301.tif 0.07 0.16 0.18 0.1810952 0.2 0.53 97553
## dimension(s):
## from to offset delta refsys point values x/y
## x 1 2919 -79.0104 0.00416667 WGS 84 FALSE NULL [x]
## y 1 2990 12.4604 -0.00416667 WGS 84 FALSE NULL [y]2.2 Atributos de un raster
## resolucion
0.00416667*110000 ## [1] 458.3337## geometria
st_bbox(luces_s)## xmin ymin xmax ymax
## -79.010415858 0.002082733 -66.847915761 12.460416166st_crs(luces_s)## Coordinate Reference System:
## User input: WGS 84
## wkt:
## GEOGCRS["WGS 84",
## DATUM["World Geodetic System 1984",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## CS[ellipsoidal,2],
## AXIS["geodetic latitude (Lat)",north,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433]],
## AXIS["geodetic longitude (Lon)",east,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433]],
## ID["EPSG",4326]]st_dimensions(luces_s)## from to offset delta refsys point values x/y
## x 1 2919 -79.0104 0.00416667 WGS 84 FALSE NULL [x]
## y 1 2990 12.4604 -0.00416667 WGS 84 FALSE NULL [y]## atributos
names(luces_s)## [1] "night_light_202301.tif"names(luces_s) = "date_202301"
## valores del raster
luces_s[[1]] %>% max(na.rm = T)## [1] 2208.49luces_s[[1]] %>% min(na.rm = T)## [1] -0.64luces_s[[1]] %>% as.vector() %>% summary() ## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## -1 0 0 0 0 2208 4031249luces_s[[1]][is.na(luces_s[[1]])==T] %>% head()# Reemplazar NA's## [1] NA NA NA NA NA NAluces_s[[1]][2000:2010,2000:2010]## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
## [1,] 0.15 0.12 0.09 0.17 0.14 0.12 0.18 0.20 0.24 0.28 0.14
## [2,] 0.15 0.12 0.07 0.16 0.14 0.13 0.12 0.16 0.16 0.25 0.14
## [3,] 0.06 0.10 0.13 0.20 0.15 0.11 0.10 0.09 0.10 0.17 0.12
## [4,] 0.07 0.13 0.18 0.16 0.11 0.20 0.21 0.08 0.12 0.14 0.12
## [5,] 0.08 0.13 0.16 0.16 0.15 0.22 0.21 0.19 0.15 0.13 0.11
## [6,] 0.12 0.13 0.16 0.15 0.13 0.18 0.19 0.17 0.15 0.10 0.13
## [7,] 0.12 0.15 0.19 0.24 0.19 0.22 0.15 0.11 0.14 0.13 0.13
## [8,] 0.16 0.11 0.08 0.18 0.21 0.19 0.11 0.15 0.13 0.19 0.19
## [9,] 0.13 0.12 0.19 0.16 0.13 0.15 0.19 0.15 0.09 0.13 0.17
## [10,] 0.15 0.15 0.16 0.16 0.11 0.06 0.07 0.14 0.11 0.11 0.18
## [11,] 0.11 0.13 0.16 0.17 0.14 0.08 0.06 0.17 0.18 0.17 0.21luces_s[[1]][2000:2010,2000:2010] %>% table() # Sustraer una parte de la matriz## .
## 0.0599999986588955 0.0700000002980232 0.0799999982118607 0.0900000035762787
## 3 3 4 3
## 0.100000001490116 0.109999999403954 0.119999997317791 0.129999995231628
## 4 10 10 15
## 0.140000000596046 0.150000005960464 0.159999996423721 0.170000001788139
## 8 14 12 7
## 0.180000007152557 0.189999997615814 0.200000002980232 0.209999993443489
## 6 9 3 4
## 0.219999998807907 0.239999994635582 0.25 0.280000001192093
## 2 2 1 1## puedo reproyectar un raster?
st_crs(luces_s)## Coordinate Reference System:
## User input: WGS 84
## wkt:
## GEOGCRS["WGS 84",
## DATUM["World Geodetic System 1984",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## CS[ellipsoidal,2],
## AXIS["geodetic latitude (Lat)",north,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433]],
## AXIS["geodetic longitude (Lon)",east,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433]],
## ID["EPSG",4326]]luces_new_crs = st_transform(luces_s,crs=4126)
luces_s[[1]][2000:2010,2000:2010] # no se alteran las geometrias## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
## [1,] 0.15 0.12 0.09 0.17 0.14 0.12 0.18 0.20 0.24 0.28 0.14
## [2,] 0.15 0.12 0.07 0.16 0.14 0.13 0.12 0.16 0.16 0.25 0.14
## [3,] 0.06 0.10 0.13 0.20 0.15 0.11 0.10 0.09 0.10 0.17 0.12
## [4,] 0.07 0.13 0.18 0.16 0.11 0.20 0.21 0.08 0.12 0.14 0.12
## [5,] 0.08 0.13 0.16 0.16 0.15 0.22 0.21 0.19 0.15 0.13 0.11
## [6,] 0.12 0.13 0.16 0.15 0.13 0.18 0.19 0.17 0.15 0.10 0.13
## [7,] 0.12 0.15 0.19 0.24 0.19 0.22 0.15 0.11 0.14 0.13 0.13
## [8,] 0.16 0.11 0.08 0.18 0.21 0.19 0.11 0.15 0.13 0.19 0.19
## [9,] 0.13 0.12 0.19 0.16 0.13 0.15 0.19 0.15 0.09 0.13 0.17
## [10,] 0.15 0.15 0.16 0.16 0.11 0.06 0.07 0.14 0.11 0.11 0.18
## [11,] 0.11 0.13 0.16 0.17 0.14 0.08 0.06 0.17 0.18 0.17 0.21luces_new_crs[[1]][2000:2010,2000:2010] # no se alteran las geometrias## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
## [1,] 0.15 0.12 0.09 0.17 0.14 0.12 0.18 0.20 0.24 0.28 0.14
## [2,] 0.15 0.12 0.07 0.16 0.14 0.13 0.12 0.16 0.16 0.25 0.14
## [3,] 0.06 0.10 0.13 0.20 0.15 0.11 0.10 0.09 0.10 0.17 0.12
## [4,] 0.07 0.13 0.18 0.16 0.11 0.20 0.21 0.08 0.12 0.14 0.12
## [5,] 0.08 0.13 0.16 0.16 0.15 0.22 0.21 0.19 0.15 0.13 0.11
## [6,] 0.12 0.13 0.16 0.15 0.13 0.18 0.19 0.17 0.15 0.10 0.13
## [7,] 0.12 0.15 0.19 0.24 0.19 0.22 0.15 0.11 0.14 0.13 0.13
## [8,] 0.16 0.11 0.08 0.18 0.21 0.19 0.11 0.15 0.13 0.19 0.19
## [9,] 0.13 0.12 0.19 0.16 0.13 0.15 0.19 0.15 0.09 0.13 0.17
## [10,] 0.15 0.15 0.16 0.16 0.11 0.06 0.07 0.14 0.11 0.11 0.18
## [11,] 0.11 0.13 0.16 0.17 0.14 0.08 0.06 0.17 0.18 0.17 0.21## plot data
plot(luces_s)## downsample set to c(7,7)
2.3 Filtrar usando la gemoetría
## download boundary
quilla <- getbb(place_name = "Barranquilla, Colombia",
featuretype = "boundary:administrative",
format_out = "sf_polygon") %>% .[1,]## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()leaflet() %>%
addTiles() %>%
addPolygons(data=quilla)## cliping
l_quilla_1 = st_crop(x = luces_s , y = quilla) # crop luces de Colombia con polygono de Medellin
l_quilla_1[[1]] %>% t()## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
## [1,] NA NA NA NA NA NA NA NA NA NA NA
## [2,] NA NA NA NA NA 32.05 72.30 NA NA NA NA
## [3,] NA NA NA NA NA NA 94.43 92.54 72.77 NA NA
## [4,] NA NA NA NA NA 95.34 76.51 56.01 59.52 47.77 NA
## [5,] NA NA NA NA NA 62.60 51.57 28.10 32.52 63.07 76.49
## [6,] NA NA NA NA 15.24 21.01 27.16 46.55 50.28 71.32 83.99
## [7,] NA NA NA NA 40.83 47.05 64.18 76.44 77.28 78.88 65.12
## [8,] NA NA NA 75.98 71.70 85.57 105.79 99.95 84.11 79.16 66.23
## [9,] NA NA NA 82.99 93.33 101.92 111.54 102.47 76.40 81.48 79.32
## [10,] NA NA NA 76.66 85.32 91.35 99.35 96.30 91.31 86.57 86.37
## [11,] NA NA 34.13 62.03 71.77 79.13 87.62 95.27 100.51 87.12 86.48
## [12,] NA NA NA 35.22 45.27 62.79 76.90 89.58 96.36 92.16 90.84
## [13,] NA NA NA 44.63 33.76 42.14 58.06 68.38 83.42 87.67 95.82
## [14,] NA NA NA 43.40 22.66 35.22 49.98 60.20 68.46 79.86 100.32
## [15,] NA NA NA NA 42.55 45.05 52.81 58.82 61.30 65.44 69.87
## [16,] NA NA NA NA 62.58 64.82 55.21 55.67 55.18 55.81 61.82
## [17,] NA 36.12 36.38 49.91 70.47 69.25 54.57 48.39 55.60 57.99 59.12
## [18,] NA 39.68 49.92 55.94 69.37 62.20 45.75 43.38 48.67 58.43 58.08
## [19,] NA NA 51.42 59.80 58.90 48.16 43.53 42.76 48.42 56.79 54.78
## [20,] NA NA 47.31 63.17 52.35 46.63 45.12 50.09 56.28 57.06 53.92
## [21,] NA NA 66.26 65.72 53.72 47.85 49.67 57.15 60.50 53.60 50.88
## [22,] NA NA 70.88 NA 56.34 52.74 53.29 67.56 65.18 53.31 55.93
## [23,] NA NA NA NA 57.27 55.93 56.22 63.81 57.03 52.76 59.40
## [24,] NA NA NA NA 54.40 56.93 57.73 57.26 51.44 49.65 50.65
## [25,] NA NA NA NA 60.24 69.57 64.83 56.72 47.53 44.60 44.62
## [26,] NA NA NA NA 54.39 69.89 67.61 56.60 48.08 45.28 47.42
## [27,] NA NA NA NA 47.97 61.82 63.79 62.72 59.50 49.69 52.51
## [28,] NA NA NA NA 25.88 43.63 49.79 63.71 60.43 52.18 59.41
## [29,] NA NA NA NA NA 47.43 54.56 61.69 58.86 60.99 69.98
## [30,] NA NA NA NA NA NA 58.19 64.35 67.04 71.71 82.92
## [31,] NA NA NA NA NA NA NA NA NA NA NA
## [,12] [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23]
## [1,] NA NA NA NA NA NA NA NA NA NA NA NA
## [2,] NA NA NA NA NA NA NA NA NA NA NA NA
## [3,] NA NA NA NA NA NA NA NA NA NA NA NA
## [4,] NA NA NA NA NA NA NA NA NA NA NA NA
## [5,] 85.56 NA NA NA NA NA NA NA NA NA NA NA
## [6,] 85.30 63.33 NA NA NA NA NA NA NA NA NA NA
## [7,] 60.87 54.69 53.39 NA NA NA NA NA NA NA NA NA
## [8,] 53.39 47.17 49.70 51.41 NA NA NA NA NA NA NA NA
## [9,] 66.51 56.38 59.51 62.07 46.97 NA NA NA NA NA NA NA
## [10,] 84.80 69.04 69.33 66.04 51.06 27.09 NA NA NA NA NA NA
## [11,] 85.51 77.50 65.75 62.99 51.63 51.20 36.54 NA NA NA NA NA
## [12,] 87.09 82.91 75.94 66.18 67.95 74.07 70.44 46.73 NA NA NA NA
## [13,] 90.66 76.26 78.83 78.96 75.55 81.83 70.24 36.41 18.83 NA NA NA
## [14,] 91.54 73.72 75.44 93.24 92.60 77.20 65.92 41.21 23.91 22.20 NA NA
## [15,] 71.60 71.90 75.94 88.63 98.86 86.89 71.64 60.35 49.21 42.73 NA NA
## [16,] 66.07 71.64 73.51 78.55 92.67 85.76 82.93 82.86 69.00 61.91 51.15 NA
## [17,] 65.07 69.80 67.73 73.24 83.59 81.23 83.38 95.06 82.41 64.15 52.33 NA
## [18,] 60.75 62.08 58.91 66.09 72.25 70.98 76.85 87.11 81.52 44.85 37.96 NA
## [19,] 54.65 60.34 60.78 67.11 64.15 64.04 64.97 68.56 65.19 51.40 72.00 NA
## [20,] 57.21 72.60 79.05 76.00 68.02 64.23 63.44 70.97 66.07 56.48 NA NA
## [21,] 54.49 66.49 70.52 66.77 71.46 68.73 68.48 75.41 76.00 55.55 NA NA
## [22,] 60.69 63.25 60.49 60.49 71.85 89.37 73.30 72.00 73.48 63.38 50.81 NA
## [23,] 61.73 75.28 80.32 70.60 82.59 84.18 70.54 71.74 76.27 80.48 83.79 NA
## [24,] 51.92 69.87 86.02 78.78 92.53 87.69 74.86 83.80 94.17 93.37 85.02 NA
## [25,] 49.16 59.19 73.68 80.36 92.06 99.48 96.92 87.09 91.57 88.35 76.77 NA
## [26,] 49.30 67.65 71.78 69.61 84.80 96.47 111.94 NA NA NA NA NA
## [27,] 58.76 80.24 76.63 67.45 72.62 79.86 97.51 NA NA NA NA NA
## [28,] 72.81 93.18 94.09 87.39 NA NA NA NA NA NA NA NA
## [29,] 79.17 96.70 97.20 NA NA NA NA NA NA NA NA NA
## [30,] 91.30 119.18 NA NA NA NA NA NA NA NA NA NA
## [31,] NA NA NA NA NA NA NA NA NA NA NA NA## Plot data
ggplot() + geom_stars(data=l_quilla_1 , aes(y=y,x=x,fill=date_202301)) + # plot raster
scale_fill_viridis(option="A" , na.value='white') +
geom_sf(data=quilla , fill=NA , col="green") + theme_bw() 
[3.] Aplicación: Actividad económica en Barranquilla
3.1 Importar datos 201301
## load data
l_quilla_0 = read_stars("input/raster/night_light_201301.tif") %>% st_crop(quilla)
names(l_quilla_0) = "date_201301"
ggplot() + geom_stars(data=l_quilla_0 , aes(y=y,x=x,fill=date_201301)) + # plot raster
scale_fill_viridis(option="A" , na.value='white') +
geom_sf(data=quilla , fill=NA , col="green") + theme_bw() 
## stack rasters: adicionar bandas
l_quilla = c(l_quilla_0,l_quilla_1)
l_quilla## stars object with 2 dimensions and 2 attributes
## attribute(s):
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## date_201301 8.70 55.10 65.90 66.55410 78.99 135.63 308
## date_202301 15.24 54.49 65.72 66.38741 78.83 119.18 308
## dimension(s):
## from to offset delta refsys point values x/y
## x 999 1021 -79.0104 0.00416667 WGS 84 FALSE NULL [x]
## y 340 370 12.4604 -0.00416667 WGS 84 FALSE NULL [y]3.2 Raster a datos vectoriales
## st_as_sf
puntos_quilla = st_as_sf(x = l_quilla, as_points = T, na.rm = T) # raster to sf (points)
puntos_quilla## Simple feature collection with 405 features and 2 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -74.85 ymin: 10.92083 xmax: -74.75833 ymax: 11.04583
## Geodetic CRS: WGS 84
## First 10 features:
## date_201301 date_202301 geometry
## 1 34.77 32.05 POINT (-74.82917 11.04167)
## 2 61.45 72.30 POINT (-74.825 11.04167)
## 3 73.49 94.43 POINT (-74.825 11.0375)
## 4 78.52 92.54 POINT (-74.82083 11.0375)
## 5 57.54 72.77 POINT (-74.81667 11.0375)
## 6 32.77 95.34 POINT (-74.82917 11.03333)
## 7 28.12 76.51 POINT (-74.825 11.03333)
## 8 27.86 56.01 POINT (-74.82083 11.03333)
## 9 36.19 59.52 POINT (-74.81667 11.03333)
## 10 45.98 47.77 POINT (-74.8125 11.03333)poly_quilla = st_as_sf(x = l_quilla, as_points = F, na.rm = T) # raster to sf (polygons)
poly_quilla## Simple feature collection with 405 features and 2 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: -74.85208 ymin: 10.91875 xmax: -74.75625 ymax: 11.04792
## Geodetic CRS: WGS 84
## First 10 features:
## date_201301 date_202301 geometry
## 1 34.77 32.05 POLYGON ((-74.83125 11.0437...
## 2 61.45 72.30 POLYGON ((-74.82708 11.0437...
## 3 73.49 94.43 POLYGON ((-74.82708 11.0395...
## 4 78.52 92.54 POLYGON ((-74.82292 11.0395...
## 5 57.54 72.77 POLYGON ((-74.81875 11.0395...
## 6 32.77 95.34 POLYGON ((-74.83125 11.0354...
## 7 28.12 76.51 POLYGON ((-74.82708 11.0354...
## 8 27.86 56.01 POLYGON ((-74.82292 11.0354...
## 9 36.19 59.52 POLYGON ((-74.81875 11.0354...
## 10 45.98 47.77 POLYGON ((-74.81458 11.0354...## Plot data
ggplot() +
geom_sf(data = puntos_quilla , aes(col=date_202301)) +
scale_color_viridis(option="A" , na.value='white') +
geom_sf(data=quilla , fill=NA , col="green") + theme_test() 
ggplot() +
geom_sf(data = poly_quilla , aes(fill=date_202301),col=NA) +
scale_fill_viridis(option="A" , na.value='white') +
geom_sf(data=quilla , fill=NA , col="green") + theme_test()
3.3 Asignar el valor del pixel a un lugar
## get polygon
puerto <- getbb(place_name = "Puerto de Barranquilla - Sociedad Portuaria, Colombia",
featuretype = "barrier:wall",
format_out = "sf_polygon")## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()leaflet() %>%
addTiles() %>%
addPolygons(data=puerto)## asignar luminosidad
l_puerto = st_join(puerto,poly_quilla)
## variacion promedio
df = l_puerto
st_geometry(df) = NULL
df %>%
summarise(pre=mean(date_201301,na.rm=T) ,
post=mean(date_202301,na.rm=T)) %>%
mutate(ratio=post/pre-1)## pre post ratio
## 1 61.6 54.46 -0.1159091[4.] Aplicación: Cambio de cobertura de suelo?
## read raster
land_cover = c(read_stars("input/raster/discreta_2015.tif"),
read_stars("input/raster/discreta_2019.tif"))
0.000992063*111000 # resolucion## [1] 110.119## read polygon
quilla_buffer <- st_buffer(x = quilla,dist = 2000)
leaflet() %>%
addTiles() %>%
addPolygons(data=quilla_buffer)quilla_border <- st_difference(x = quilla_buffer , y = quilla)
leaflet() %>%
addTiles() %>%
addPolygons(data=quilla_border)## cliping raster
quilla_cover = land_cover %>% st_crop(quilla_buffer)
plot(quilla_cover) ## view data
## cliping raster
border = land_cover %>% st_crop(quilla_border)
plot(border) ## view data
## raster to sf
border_sf = st_as_sf(x=border, as_points = F, na.rm = T)
## cuantos pixeles se urbanizaron?
border_sf = border_sf %>%
rename(c_2015=discreta_2015.tif,c_2019=discreta_2019.tif)
border_sf = border_sf %>%
mutate(build = ifelse(c_2015!=50 & c_2019==50,1,0))
table(border_sf$build)##
## 0 1
## 7970 52[5.] Aplicación: Clusteres espaciales
## restaurantes
points <- opq(bbox = getbb("Bogota Colombia")) %>%
add_osm_feature(key = "amenity", value = "restaurant") %>%
osmdata_sf() %>% .$osm_points %>% select(osm_id,name) %>% mutate(restaurant=1)
leaflet() %>% addTiles() %>% addCircles(data=points)## download boundary
city <- getbb(place_name = "Bogota, Colombia",
featuretype = "boundary:administrative",
format_out = "sf_polygon") %>% .[[2]]## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()leaflet() %>% addTiles() %>% addPolygons(data=city)## afine transformations
points = st_transform(points,3857)
city = st_transform(city,3857)
## sf to ppp
points_p <- as.ppp(X = points["restaurant"])
Window(points_p) <- as.owin(city)
plot(points_p)
#--------------#
## summary
summary(points_p) # las unidades son en metros## Marked planar point pattern: 6751 points
## Average intensity 1.697305e-05 points per square unit
##
## Coordinates are given to 1 decimal place
## i.e. rounded to the nearest multiple of 0.1 units
##
## marks are numeric, of type 'double'
## Summary:
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1 1 1 1 1 1
##
## Window: polygonal boundary
## 5 separate polygons (no holes)
## vertices area relative.area
## polygon 1 7588 395800000.0 0.995000
## polygon 2 276 1524690.0 0.003830
## polygon 3 72 235244.0 0.000591
## polygon 4 45 75844.4 0.000191
## polygon 5 81 112582.0 0.000283
## enclosing rectangle: [-8262524, -8238784] x [498234.9, 538665.2] units
## (23740 x 40430 units)
## Window area = 397748000 square units
## Fraction of frame area: 0.414intensity(points_p)## [1] 1.697305e-05## reescalar
points_p = rescale(X = points_p , s = 100 , unitname = "100-metros")
unitname(points_p)## 100-metros / 100-metrossummary(points_p)## Marked planar point pattern: 6751 points
## Average intensity 0.1697305 points per square 100-metros
##
## Coordinates are given to 3 decimal places
## i.e. rounded to the nearest multiple of 0.001 100-metros
##
## marks are numeric, of type 'double'
## Summary:
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1 1 1 1 1 1
##
## Window: polygonal boundary
## 5 separate polygons (no holes)
## vertices area relative.area
## polygon 1 7588 39580.00000 0.995000
## polygon 2 276 152.46900 0.003830
## polygon 3 72 23.52440 0.000591
## polygon 4 45 7.58444 0.000191
## polygon 5 81 11.25820 0.000283
## enclosing rectangle: [-82625.24, -82387.84] x [4982.349, 5386.652] 100-metros
## (237.4 x 404.3 100-metros)
## Window area = 39774.8 square 100-metros
## Unit of length: 1 100-metros
## Fraction of frame area: 0.414intensity(points_p)## [1] 0.1697305## estimated standard error for intensity
sqrt(intensity(points_p)/summary(points_p)$window$area)## [1] 0.002065741## plot
plot(density(points_p,5)) ## plot density
contour(density(points_p,bw.ppl(points_p)), axes = FALSE) ## contour density
#--------------#
## quadrants
points_q = quadratcount(points_p, nx = 4, ny = 6)
points_q## tile
## Tile row 1, col 3 Tile row 1, col 4 Tile row 2, col 2 Tile row 2, col 3
## 0 20 110 335
## Tile row 2, col 4 Tile row 3, col 1 Tile row 3, col 2 Tile row 3, col 3
## 759 0 362 835
## Tile row 3, col 4 Tile row 4, col 1 Tile row 4, col 2 Tile row 4, col 3
## 806 54 650 2031
## Tile row 4, col 4 Tile row 5, col 1 Tile row 5, col 2 Tile row 5, col 3
## 237 13 187 341
## Tile row 5, col 4 Tile row 6, col 2 Tile row 6, col 3
## 0 0 11plot(points_q, cex = 2 , col="red")
Referencias
Lovelace, R., Nowosad, J., & Muenchow, J. (2019). Geocomputation with R. [Ver aquí]
- Cap. 4: Spatial Data Operations
- Cap. 5: Geometry Operations
- Cap. 6: Reprojecting geographic data
- Cap. 11: Statistical learning
Bivand, R. S., Pebesma, E. J., Gómez-Rubio, V., & Pebesma, E. J. (2013). Applied spatial data analysis with R. [Ver aquí]
- Cap. 4: Spatial Data Import and Export
- Cap. 7: Spatial Point Pattern Analysis
- Cap. 8: Interpolation and Geostatistics